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Abstract. The boundary integral method (BIM) is a formulation of Helmholtz 
equation in the form of an integral equation suitable for numerical discretiza- 
tion to solve the quantum billiard. This paper is an extensive numerical sur- 
vey of BIM in a variety of quantum billiards, integrable (circle, rectangle), 
KAM systems (Robnik billiard) and fully chaotic (ergodic, such as stadium, 
Sinai billiard and cardioid billiard). On the theoretical side we point out some 
serious flaws in the derivation of BIM in the literature and show how the final 
formula (which nevertheless was correct) should be derived in a sound way 
and we also argue that a simple minded application of BIM in nonconvex 
geometries presents serious difficulties or even fails. On the numerical side 
we have analyzed the scaling of the averaged absolute value of the systematic 
error AE of the eigenenergy in units of mean level spacing with the density 
of discretization (b = number of numerical nodes on the boundary within one 
de Broglie wavelength), and we find that in all cases the error obeys a power 
law < \AE\ >= Ab~ a , where a (and also A) varies from case to case (it is 
not universal), and is affected strongly by the existence of exterior chords in 
nonconvex geometries, whereas the degree of the classical chaos seems to be 
practically irrelevant. We comment on the semiclassical limit of BIM and 
make suggestions about a proper formulation with correct semiclassical limit 
in nonconvex geometries. 
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1 Introduction 



In the studies of quantum chaos (Gutzwiller 1990, Giannoni et al 1989) good 
numerical methods are not only indispensable but quite essential. We need 
them in order to illustrate and verify the theoretical developments. Also, new 
numerical experiments provide valuable material as evidence and inspiration 
for new theories. Quantum billiards are certainly very useful model systems 
since dynamically they are generic and - depending on the design - they can 
cover all regimes of classical motion between integrability and full chaos (er- 
godicity). Their classical dynamics can be easily followed for very long time 
periods because no integration is necessary but only searching for zeros of 
certain functions at collision points. Quantanlly they have the advantage 
of having a compact configuration space admitting application of various 
numerical methods such as the Boundary Integral Method (BIM) (Banerjee 
1994, Berry and Wilkinson 1984, henceforth referred to as BW, and the refer- 
ences therein), the Plane Wave Decomposition Method (PWDM) employed 
by Heller (1984), conformal mapping diagonalization technique introduced 
by Robnik (1984) and further developed by Berry and Robnik (1986) and 
recently by Prosen and Robnik (1993, 1994), also by (Bohigas et al 1993), 
and a number of other methods, which, however, might be adapted to some 
special systems. Among the general methods BIM is probably most widely 
used, even in quite practical engineering problems. The main task of the 
present paper is to point out some logical flaws in deriving this method 
in the context of the existing literature, to analyze its limitations in cases 
of nonconvex geometries with suggestions for improvements and generaliza- 
tions, and to perform an extensive numerical investigation of its accuracy in 
relation to geometrical properties and classical dynamics. 

In order to clearly expose the difficulties and the errors in the derivation 
of BIM offered in the literature, e.g. in BW, we present our regularized 
derivation, by which we mean that we construct and use a Green function 
which automatically (by construction) satisfies the Dirichlet boundary con- 
dition (vanishes on the boundary dB of the billiard domain £>), which is 
achieved by employing the method of images (see e.g. Balian and Bloch 
(1974) and the references therein). This will enable us to avoid committing 
two errors, which, however, luckily compensated each other: Firstly, in tak- 
ing the normal derivatives on the two sides of equation (6) in BW, on the rhs 
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we must use the value of i[)(r) which is the interior solution inside B rather 
than |V( r ) which is the value exactly on the boundary, simply because in 
taking the derivatives we must evaluate the function at two infinitesimally 
separated points normal to the boundary; Secondly, this error of taking the 
unjustified factor 1/2 is then exactly compensated by another error in arriv- 
ing at the equation (8) in BW, namely by interchanging the integration along 
the boundary dB and the normal differentiation, because due to singularities 
on the boundary these two operations do not commute. 

So, let us offer our regularized derivation. We are searching for the solu- 
tion ^>(r) with eigenenergy E = k 2 obeying the Helmholtz equation 

V^(r) + ^(r) = 0, (1) 

with the Dirichlet boundary condition ^>(r) =0 on the boundary r G dB. We 
will transform this Schrodinger equation for our quantum billiard B into an 
integral equation by means of the regularized Green function G(r, r'), which 
solves the following defining equation 

V*G(r, r') + k 2 G{r, r') = 5{r - r') - 6{r - r' R ), (2) 

where r and r' are in B U dB, and t'r is the mirror image of r' with respect 
to the tangent at the closest lying point on the boundary, and thus if r' is 
sufficiently close to the boundary then t'r is outside the billiard B. The 
solution can easily be found in terms of the free propagator (the free particle 
Green function on the full Euclidean plane) 

G (r,r') = -^ 1) (A;|r-r'|), (3) 

where Hq is the zero order Hankel function of the first kind (Abramowitz 
and Stegun 1972), namely 

G(r,r / ) = Go(r,r / )-G (r,r / i? ), (4) 

such that now G(r, r') is zero by construction for any r' on the boundary 
dB, in contradistinction to the Green function defined and used in equation 
(5) in BW. Multiplication of the equation (0) by ip(r) and the Helmholtz 
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equation ([I]) by G(r,r'), subtraction, integration over the area inside B and 
using Green's theorem, yields 



j ds (V>(r)n • V r G(r, r') - G(r, r')n • V r ^(r)) = i/;(r'), (5) 

where s is the arclength on the boundary dB oriented anticlockwise, n is 
the unit normal vector to dB at r oriented outward, and this equation is 
now valid for all r' inside and on the boundary of B. Since in this equation 
everything is regular we can take the normal partial derivatives on both sides. 
Following the usual notation in BW we define the normal derivative of i/j at 
the point s as 

u(s) = n- V r V>(r(s)), (6) 
and thus using the boundary condition ip(r) = we arrive at 

u(s) = -2 j ds'u(s')n ■ V r G (r, r'), (7) 

In this way we have correctly derived the main integral equation of the bound- 
ary integral method which is correctly given as equation (8) in BW (where 
the two errors exactly compensate), so that all the further steps in working 
out the geometry of equation (0) and the numerical discretization are exactly 
the same as in BW. As shown in figure 1 we define the length of the chord 
between two points on the boundary r(s) and r(s') as 

p(s,s') = \r(s)-r(s')\ (8) 

and the angle 6(s', s) is the angle between the chord and the tangent to dB 
at s'. Of course 6(s, s') ^ @{s', s). Thus following the notation of BW we can 
write 

and we obtain finally 

u{s) = - l -ik j ds'u(s') sin#(s, s^H^ {kp(s, s')}. (10) 

In numerically solving this integral equation we have used precisely the same 
primative discretization procedure as BW, which turned out to be better than 
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some other more sophisticated versions. So we simply divide the perimeter 
C into iV equally long segments and thus define 



s m = mC/N, p(s t , s m ) = pim, 9(si, s m ) = 9 lm , l<l,m<N, (11) 

Therefore numerically we are searching for the zero of the determinant A N (E) = 
det(Mi m ) where M[ m are the matrix elements of the N x N matrix, 



where E = k 2 . Due to the asymmetry 0i m ^ m i this matrix is a general 
complex non-Hermitian matrix. 

One important aspect of this formalism is the semiclassical limiting form 
which has been extensively studied by Boasman (1994) and which is the sub- 
ject of our current investigation (Li and Robnik 1995a). At this point we 
just want to make the following comment. In cases of nonconvex geome- 
tries we will have exterior chords connecting two points on the boundary 
such that they lie entirely or at least partially, outside B. While formally 
the method and the procedure in such cases is perfectly right, in reality it 
fails completely, and this can be seen by considering the semiclassical limit. 
The formal leading order in the asymtotic expansion of the Hankel function 
(Debye approximation) does not match the actually correct semiclassi- 
cal leading approximation which would be spanned by the shortest classical 
orbit connecting the two points via at least one or many collision points. 
Therefore we understand and expect that the method must fail or at least 
must meet severe difficulties in cases of nonconvex geometry. This has been 
completely confirmed in our present work as we will show in the next section, 
whilst the analytical work to reformulate the method including the multiple 
collision expansions is in progress (Li and Robnik 1995b) and is expected to 
deal satisfactorily with nonconvex geometries. 

In the next section we shall analyze the numerical accuracy of BIM as a 
function of the density of discretization 



Mi m = Sim + 77T7 sin Ol m H{ \kpi m ), 



(12) 



b = 



2nN 



(13) 



kC ' 
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in a variety of quantum billiards with integrable, KAM-type or ergodic clas- 
sical dynamics, including such with nonconvex geometry. The main result 
is that there is always a power law so that the error of eigenenergy in units 
of mean level spacing, after taking average of the absolute value over a suit- 
able ensemble of eigenstates, obeys < \AE\ >= Ab~ a , but the exponent a is 
nonuniversal and typically becomes almost zero if there are nonconvex seg- 
ments on the boundary. 



2 Numerical results 

The numerical procedure we have used to solve the variety of quantum bil- 
liards is exactly as described in the introduction and therefore it is precisely 
the same as in BW. Our main task in this work is to analyze in detail the 
behavior of BIM as a function of the density of discretization b, especially in 
relation to the geometrical properties of B (nonconvexities) and in relation to 
classical dynamics, whose chaotic behaviour is expected to imply interesting 
methodological and algorithmical manifestation of quantum chaos. 

This to end we have to measure the numerical error of the eigenenergies 
in some natural units, which obviously is the mean level spacing. In plane 
billiards this is well defined and determined by the leading term of the Weyl 
formula, namely in our units it is equal to 4ir/A, where A is the area of the 
billiard B. Since the error AE of the eigenenergies thus defined still fluctu- 
ates widely from state to state we have to perform some kind of averaging 
over a suitable ensemble of consecutive states. But this will make sense only 
if such a local average of the error is stationary (constant) over a suitable 
energy interval. This condition has been confirmed to be satisfied in almost 
all cases which we checked. In case of the circle billiard and in case of half 
circle billiard (which embodies all the odd states of the full circle billiard) 
this stationarity of the locally averaged error is shown in figures 2(a-b), re- 
spectively. 

Having established that we have always taken the average of the absolute 
value of the error (in units of the mean level spacing) over a suitable ensemble 
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of eigenstates, for which we have chosen the lowest one hundred eigenstates 
in all cases. This quantity will be henceforth denoted by < |Ai?| >. 

In figure 3(a-f) we show the error < \AE\ > versus b for three shapes of 
the Robnik billiard with the shape parameter A = 0, 1/4, 1/2. (The billiard 
shape is defined by the quadratic conformal map of the unit disk \z\ < 1 onto 
to the w-plane, namely w(z) = z + \z 2 .) We show the normal plot of the 
averaged absolute value of the error versus b in the figures 3(a,c,e) and the 
log- log plot in figures 3(b,d,f). The best fitting power law curve is described 
by 

< |A£7| >= Ab~ a , (14) 

and is seen to provide a very significant fit in all three cases. One should 
be reminded that at A = we have integrable classical dynamics in the 
circle billiard, at A = 1/4 we have almost ergodic but nevertheless KAM- 
type dynamics with very tiny islands of stability (see e.g. Li and Robnik 
(1994,1995c); it should be emphasized that at A = 1/4 we have zero curva- 
ture point at z = —1, and for all A > 1/4 the shape is nonconvex), whilst 
at A = 1/2 we have rigorous ergodicity (Markarian 1993) and also noncon- 
vex geometry. So one can observe that increasing chaos from integrability 
(circle) to almost ergodicity (A = 1/4) has almost no effect on a, whereas 
the nonconvexities at A > 1/4 seem to imply a complete "collapse" of a: 
As it will be shown in figure 5(a) a drops to zero almost discontinuously at 
A = 1/4 and then increases up to 0.4 at A close to 1/2. This is because, 
perhaps, close to A = 1/2 where we have the cusp singularity at z = — 1, 
the role of the nonconvexities is smaller. In all cases we have calculated all 
states by applying BIM but then for technical reasons compared only the odd 
states with their exact value, which are supplied by the conformal mapping 
diagonalization technique (see e.g. Robnik 1984, Prosen and Robnik 1993). 

It is then interesting to look at the similar plots for BIM as applied to half 
billiard by which we mean the upper part of the Robnik billiard defined for 
^s(z) > 0, implying also Q(w) > 0, still with the Dirichlet boundary condi- 
tions everywhere on the boundary. This embodies all the odd states of the 
full billiard. The classical dynamics in the two billiards is of course exactly 
the same and yet in figures 4(a-f) for the half billiard we see that a is now 
notably different, showing that there is no universality in the value of a so 
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that a is certainly not uniquely determined by the classical dynamics of the 
underlying quantum billiard. On the other hand we can also see the effect of 
nonconvexities of the boundary OB: Namely for the half (cardioid) billiard at 
A = 1/2 we have no nonconvexities and this probably is precisely the reason 
for large increase of a from 0.4 in figure 3(e-f) to 2.4 in figure 4(e-f). 

Having established the validity of the power law flT4]) it is now most interest- 
ing and also immensely CPU-time consuming (it took almost one month of 
CPU-time on Convex C3860 to produce figure 5(a,c) and a little bit less for 
figure 5(b,d)) to look at the variation of a with the billiard shape parameter 
A. For the full Robnik billiard this is shown in figure 5(a) and for the half 
billiard in figure 5(b). In both cases there is a flat region of almost constant 
a within < A < 1/4: In the former case it fluctuates slightly around 3.5 
and in the latter case it is surprisingly stable around 2.9. At A > 1/4 the 
nonconvexities of boundary appear and this implies - as explained qualita- 
tively in the introduction - strong drop of a which is much sharper (almost 
discontinuous) in case of the full billiard because obviously the nonconvex- 
ities are more pronounced there than in the half billiard. So in the former 
case a drops almost to zero, whereas in the half billiard it decreases rather 
smoothly to about 0.16 at A = 0.36, and then increases again reaching the 
value of 2.4 at A = 1/2. 

Apart from a in (|14|) we would also like to know the value of the constant A 
(the pre-factor) in each case. This is given by fixing b = 12 and plotting the 
mean absolute value of the error (averaged over the lowest one hundred odd 
states) versus A logarithmically. Here we see that in both cases the mean er- 
ror < \ AE\ > is almost constant up to A < 1/4 and is about 5 x 10~ 6 for the 
full billiard whilst for the half billiard it is about 7 x 10 -5 . This discrepancy 
by almost an order of magnitude is not completely understood. At A > 1/4 
please observe the dramatic increasing of < |Ai£| >. For the full billiard it 
increases almost discontinuously up to 2 x 10~ 2 whereas in the half billiard 
it reaches the minimum of 4 x 10~ 3 at about A = 0.36 and decreases then 
again to 2 x 10~ 4 at A = 1/2. This difference again is qualitatively explained 
by the role of nonconvexities. 

Finally in figure 6(a-f) we show the worst cases of applying BIM, namely 
the Robnik billiard at A = 0.4 in figures 6(a,b), the half billiard in figures 
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6(c,d) with the same value of A = 0.4 and the desymmetrized Sinai billiard 
(1/8 of the full billiard) with the unit side length and radius R = 1/4 of the 
circular obstacle. In all cases a is almost zero or even slightly negative which 
implies that there is practically no convergence of the numerical eigenvalues 
by increasing the density of discretization b. Clearly, this is due to the non- 
convex geometry of these billiards. 

Most of our results are summarized in table 1 for three classes of billiard 
systems with different type of classical dynamics, namely integrable, KAM- 
type and ergodic systems. We show the calculated values of a and also the 
average absolute value of the error < \AE\ > with fixed value of b = 12. 
The table clearly demonstrates that the power law ([TJ]) for BIM is univer- 
sal but not the exponent a and the prefactor A. It also demonstrates that 
classical dynamics has little effect on a whereas the nonconvexities of the 
boundary are quite crucial: In all cases of pronounced nonconvex geometry 
a is typically close to zero or even slightly negative like in the Sinai billiard. 
In the table we include also the results on the integrable case of the rectan- 
gular equilateral triangle (half of the unit square) where a = 3.28 is quite 
large, and the ergodic case of the 1/4 Heller's stadium (2x2 square plus 
two semicircles with unit radius) in which case a = 3.0 is also quite large. 
Both of these two cases have convex geometry and thus large a despite the 
completely different classical dynamics. 

When thinking about improving the efficiency and the accuracy of BIM we 
have also tried a more sophisticated version of BIM, where we have explic- 
itly used a Gaussian integration on the boundary when discretizing our main 
equation fl7|). However, this experience has been negative after many carefull 
checks in various billiards, and therefore we decided to resort to the primitive 
discretization of BIM which is exactly the same approach as in BW. 

3 Discussion and conclusions 

The main purpose of this paper is to investigate the behaviour of the Bound- 
ary Integral Method (BIM) with respect to the density of discretization b as 
defined in equation ( |13"D (b = the number of numerical nodes per de Broglie 
wavelength along the boundary). In all cases we discovered that there is a 
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power law behaviour described in (0): The average absolute value of the er- 
ror measured in the units of mean level spacing is given by < \AE\ >= Ab~ a . 
We wanted to verify whether there is any systematic effect of classical dynam- 
ics of quantum billiards on a and A. The answer is negative. On the other 
hand we found that the role of nonconvex geometry of the boundary is cru- 
cial: If there is a nonconvex part of the boundary a is typically close to zero 
or even slightly negative implying that there is practically no convergence of 
the eigenvalues with respect to the increasing b. This failure of BIM can be 
understood as explained in the introduction, and the easiest way to see that 
is to consider the semiclassical limiting approximation of BIM (Li and Rob- 
nik 1995a). After explaining two systematic errors in the literature where the 
integral BIM equation is derived and where luckily the two errors mutually 
exactly compensate, we have given the correct (regularized) derivation and 
discussed the BIM formalism thus derived. We agree that even in nonconvex 
geometries it is formally right, but nevertheless practically fails which is most 
clearly demonstrated by the semiclassical limit mentioned above. Therefore 
we suggest a generalization of the BIM method by using multiple reflection 
(collision) expansion which is another subject of our current investigation (Li 
and Robnik 1995b), and is important not only for studies in quantum chaos 
but also in engineering problems (Banerjee 1994). 

Most of our results are summarized in table 1, giving the evidence for the 
above conclusions. In another work (Li and Robnik 1995d) we have investi- 
gated the impact of classical chaos on another general numerical method for 
quantum billiards, namely the plane wave decomposition method employed 
by Heller (1984, 1991), where we have found that at fixed b the average ab- 
solute value of the error < \AE\ > does correlate with classical chaos and 
increases sharply with increasing classical chaos. In this method such a be- 
haviour can be understood much more easily: In classically integrable cases 
the wavefunction in the semiclassical limit can be correctly described locally 
by a finite number of plane waves, whereas in classically fully chaotic (er- 
godic) systems we need locally infinite number of plane waves. Therefore at 
fixed b the accuracy of Heller's method strongly deteriorates as the system 
approaches ergodicity. More results on that will be published in a separate 
paper (Li and Robnik 1995d). 

It remains an interesting and important theoretical problem to study the 
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sensitivity of the eigenstates (eigenenergies and wavefunctions) on the bound- 
ary data of eigenfunctions, of which one aspect is also the dependence of the 
eigenstates on the billiard shape parameter. If such sensitivity correlates with 
classical chaotic dynamics and at the same time manifests itself in the accu- 
racy of the purely quantal numerical methods then such a behaviour would 
be one important manifestation of quantum chaos. This interesting line of 
thoughts in the search of another face of quantum chaos will be further devel- 
oped in another work (Li and Robnik 1995e) where we also present detailed 
studies of level curvature distribution and other measures of the sensitivity 
of the eigenstates. 
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Table 



Table 1. The power law exponent a and the average absolute value of the 
error < \AE\ > with b = 12 for different billiards. For details of the KAM- 
type see also figures 5(a,b). 



Type 


Quantum billiard 


a 


< \AE\ > 6 =12 




Circle (half) 


2.94 ± 0.17 


6.74E-5 


Integrable 


Circle (full) 


3.44 ± 0.18 


5.97E-6 




Rectangle-triangle 


3.28 ± 0.29 


4.08E-5 


KAM 


Robnik (full) 

(0 < A < 1/4) 


« 3.4 


« 5.0E-6 




Robnik (half) 

(0 < A < 1/4) 


« 2.9 


« 7.0E-5 




Stadium (1/4) 


3.00 ± 0.16 


1.18E-4 




Cardioid (half) 


2.42 ± 0.11 


1.76E-4 




Cardioid (full) 


0.42 ± 0.08 


1.85E-2 


Ergodic 


Sinai (1/8) 


- 0.34 ± 0.11 


3.63E-1 




Robnik (full) 

(0.3 < A < 1/2) 


see figure 5 a 


see figure 5 c 




Robnik (half) 
(0.3 < A < 1/2) 


see figure 5b 


see figure 5d 
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Figure captions 



Figure 1: The notation of the angles and chords used in BIM. 

Figure 2: The BIM error (measured in units of mean level spacing) of 
eigenstates versus energy. The error is difference between the BIM value and 
the exact value. The lowest 1000 odd states are shown. Plot (a) is for the 
full circle billiard and (b) for the half circle billiard. In both cases b is fixed, 
6 = 6. 

Figure 3: The ensemble averaged (over 100 lowest odd eigenstates) absolute 
BIM error versus the density of boundary discretization b and the best power 
law fit for the full Robnik billiard at different shape parameters A. The + 
represents the numerical data, the curve is the best power law fit whose a is 
given in (b,d,f). Figures (b,d,f) are log-log plots of the same quantaties as in 
(a,c,e). 

Figure 4: The same as in figure 3 but for the half Robnik billiard. 

Figure 5: We show a versus A plot for Robnik billiard in (a) and for half 
Robnik billiard in (b). In (c) and (d) we plot — lg(< |Ai?| >) with fixed 
b = 12 versus A for the two billiards, respectively. 

Figure 6: We show three examples of worst cases in applying BIM, Robnik 
billiard with A = 0.4 in (a,b), half Robnik billiard with the same A in (c,d), 
and desymmetrized Sinai billiard (1/8) with unit side length and radius R = 
1/4. In (a,c,e) we plot the averaged absolute value of the error versus b, and 
in (b,d,f) we plot log-log diagrams of the same quantities. 
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